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Abstract 

We propose a method for the construction of preconditioners of parameter-dependent 
matrices for the solution of large systems of parameter-dependent equations. The pro¬ 
posed method is an interpolation of the matrix inverse based on a projection of the 
identity matrix with respect to the Frobenius norm. Approximations of the Frobenius 
norm using random matrices are introduced in order to handle large matrices. The re¬ 
sulting statistical estimators of the Frobenius norm yield quasi-optimal projections that 
are controlled with high probability. Strategies for the adaptive selection of interpola¬ 
tion points are then proposed for different objectives in the context of projection-based 
model order reduction methods: the improvement of residual-based error estimators, 
the improvement of the projection on a given reduced approximation space, or the 
re-use of computations for sampling based model order reduction methods. 


1 Introduction 

This paper is concerned with the solntion of large systems of parameter-dependent eqnations 
of the form 

n(0«(0 = KO, (1) 

where ^ takes valnes in some parameter set S. Snch problems occnr in several contexts 
snch as parametric analyses, optimization, control or uncertainty quantihcation, where ^ 
are random variables that parametrize model or data uncertainties. The efficient solution 
of equation Q generally requires the construction of preconditioners for the operator A(^), 

iEcole Centrale de Nantes, GeM, UMR CNRS 6183, France. 

Corresponding author (anthony.nouy@ec-nantes.fr). 

*This work was supported by the French National Research Agency (Grant ANR CHORUS MONU-0005) 


1 



either for improving the performance of iterative solvers or for improving the quality of 
residual-based projection methods. 

A basic preconditioner can be dehned as the inverse (or any preconditioner) of the matrix 
A(^) at some nominal parameter value ^ G S or as the inverse (or any preconditioner) of a 
mean value of A(,^) over S (see e.g. [211 ED])- When the operator only slightly varies over 
the parameter set S, these parameter-independent preconditioners behave relatively well. 
However, for large variabilities, they are not able to provide a good preconditioning over 
the whole parameter set S. A hrst attempt to construct a parameter-dependent precondi¬ 
tioner can be found in HZ!, where the authors compute through quadrature a polynomial 
expansion of the parameter-dependent factors of a LU factorization of A(,^). More recently, 
a linear Lagrangian interpolation of the matrix inverse has been proposed in m- The gen¬ 
eralization to any standard multivariate interpolation method is straightforward. However, 
standard approximation or interpolation methods require the evaluation of matrix inverses 
(or factorizations) for many instances of ^ on a prescribed structured grid (quadrature or 
interpolation), that becomes prohibitive for large matrices and high dimensional parametric 
problems. 

In this paper, we propose an interpolation method for the inverse of matrix A(^). The 
interpolation is obtained by a projection of the inverse matrix on a linear span of samples of 
A(^)“^ and takes the form 

m 

i=l 

where are m arbitrary interpolation points in S. A natural interpolation could 

be obtained by minimizing the condition number of Pm(0^(0 -^*(0) which is a 

Clarke regular strongly pseudoconvex optimization problem [30] • However, the solution of 
this non standard optimization problem for many instances of ^ is intractable and proposing 
an efficient solution method in a multi-query context remains a challenging issue. Here, 
the projection is dehned as the minimizer of the Frobenius norm of / — that 

is a quadratic optimization problem. Approximations of the Frobenius norm using random 
matrices are introduced in order to handle large matrices. These statistical estimations of 
the Frobenius norm allow to obtain quasi-optimal projections that are controlled with high 
probability. Since we are interested in large matrices, A(^i)“^ are here considered as implicit 
matrices for which only efficient matrix-vector multiplications are available. Typically, a 
factorization (e.g. LU) of A(^j) is computed and stored. Note that when the storage of 
factorizations of several samples of the operator is unaffordable or when efficient precondi¬ 
tioners are readily available, one could similarly consider projections of the inverse operator 
on the linear span of preconditioners of samples of the operator. However, the resulting 
parameter-dependent preconditioner would be no more an interpolation of preconditioners. 
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This straightforward extension of the proposed method is not analyzed in the present paper. 

The paper then presents several contribntions in the context of projection-based model 
order redaction methods (e.g. Rednced Basis, Proper Orthogonal Decomposition (POD), 
Proper Generalized Decompositon) that rely on the projection of the solntion u{^) of ([^ on 
a low-dimensional approximation space. We hrst show how the proposed preconditioner can 
be nsed to define a Galerkin projection-based on the preconditioned residnal, which can be 
interpreted as a Petrov-Galerkin projection of the solntion with a parameter-dependent test 
space. Then, we propose adaptive constrnction of the preconditioner, based on an adap¬ 
tive selection of interpolation points, for different objectives: (i) the improvement of error 
estimators based on preconditioned residnals, (ii) the improvement of the qnality of projec¬ 
tions on a given low-dimensional approximation space, or (iii) the re-nse of compntations for 
sample-based model order redaction methods. Starting from a m-point interpolation, these 
adaptive strategies consist in choosing a new interpolation point based on different criteria. 
In (i), the new point is selected for minimizing the distance between the identity and the 
preconditioned operator. In (ii), it is selected for improving the qnasi-optimality constant of 
Petrov-Galerkin projections which measnres how far the projection is from the best approx¬ 
imation on the rednced approximation space. In (iii), the new interpolation point is selected 
as a new sample determined for the approximation of the solntion and not of the operator. 
The interest of the latter approach is that when direct solvers are nsed to solve equation ([^ 
at some sample points, the corresponding factorizations of the matrix can be stored and the 
preconditioner can be computed with a negligible additional cost. 

The paper is organized as follows. In Section]^ we present the method for the interpolation 
of the inverse of a parameter-dependent matrix. In Section]^ we show how the preconditioner 
can be used for the definition of a Petrov-Galerkin projection of the solution of ([^ on a given 
reduced approximation space, and we provide an analysis of the quasi-optimality constant 
of this projection. Then, different strategies for the selection of interpolation points for the 
preconditioner are proposed in Section]^ Finally, in Section]^ numerical experiments will 
illustrate the efficiency of the proposed preconditioning strategies for different projection- 
based model order reduction methods. 

Note that the proposed preconditioner could be also used (a) for improving the quality 
of Galerkin projection methods where a projection of the solution u{^) is searched on a 
subspace of functions of the parameters (e.g. polynomial or piecewise polynomial spaces) 
[IS1EI1E3], or (b) for preconditioning iterative solvers for ([^, in particular solvers based on 
low-rank truncations that require a low-rank structure of the preconditioner [2H1 |32l |22l [23] . 
These two potential applications are not considered here. 
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2 Interpolation of the inverse of a parameter-dependent 
matrix using Probenius norm projection 

In this section, we propose a construction of an interpolation of the matrix-valued function 
^ !-)■ G for given interpolation points in 2. We let P* = 

1 < i < m. For large matrices, the explicit computation of Pj is usually not affordable. 
Therefore, P* is here considered as an implicit matrix and we assume that the product of P* 
with a vector can be computed efficiently. In practice, factorizations of matrices A('Ci) are 
stored. 


2.1 Projection using Probenius norm 

We introduce the subspace Ym = spanjPi,... ,Pm} of An approximation Pm(0 of 

in Ym is then dehned by 

PmiO = argmin ||J - PA(0 ||f, (2) 

P&Ym 

where I denotes the identity matrix of size n, and || ■ \\f is the Probenius norm such that 
||P|||^ = (P, B)f with (P, C)f = trace(P^C'). Since G Wi, we have the interpolation 

property Pm{^i) = 1 < i < m. The minimization of ||/ — PA||i7’ has been hrst 

proposed in [25] for the construction of a preconditioner P in a subspace of matrices with 
given sparsity pattern (SPAI method). The following proposition gives some properties of 
the operator Lemma 2.6 and Theorem 3.2 in [M])- 

Proposition 2.1 Let Pm{0 defined by (|^. IFe have 

(1 - amior <\\I- Pn.mml < n{l - aim (3) 


where am(0 lowest singular value of PmiO^iO verifying 0 < am{0 — 3? 

Pm(0^(0 = I if and only if am{D = 1- Also, the following bound holds for the condi¬ 
tion number o/Pm(0^(0-' 


KiPrniOm) < 


Vn - {n - l)al{^) 

amiO 


Under the condition ||/ — Pm(0^(0lli^ < 3; eguations ([^ and Q imply that 


n{Pm{0A{0) < 


^n-{n- 1)(1 - ||/ - P^(0A(0I;F 

l-||/-P^(e)A(0||F 


( 4 ) 
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For all A G M"*, we have 


IIJ - X^P^Aml = n- 2X'^S{0 + X^MiOX, 

i=l 

where the matrix M(^) G and the vector S'(^) G are given by 

= trace(A'^(0^^^j^(0) and S'i(0 = trace(Piv4(^)). 

Therefore, the solntion of problem (|^ is PmiO = -^*(0-^* with A(^) the solntion of 

M(^)A(^) = *S'(^). When considering a small nnmber m of interpolation points, the compn- 
tation time for solving this system of eqnations is negligible. However, the compntation of 
M(^) and S'(^) reqnires the evalnation of traces of matrices AA{^)Pf PjA{^) and PiA{^) for 
all 1 < i,j < m. Since the P* are implicit matrices, the compntation of snch prodncts of 
matrices is not affordable for large matrices. Of conrse, since trace(P) = 
trace of an implicit matrix B conld be obtained by compnting the prodnct of B with the 
canonical vectors ei,..., e„, bnt this approach is clearly not affordable for large n. 

Hereafter, we propose an approximation of the above constrnction nsing an approximation 
of the Frobenins norm which reqnires less compntational efforts. 

2.2 Projection using a Frobenins semi-norm 

Here, we define an approximation Pm{0 of by 

PUO = argmin ||(/ - PA{0)V\\f, (5) 

P&Ym 

where V G with K < n. P i—>■ ||PH||i;’ dehnes a semi-norm on Here, we assnme 

that the linear map P PA{^)V is injective on so that the solntion of ([^ is nniqne. 
This reqnires K > m and is satished when rank(l/) > m and is the linear span of linearly 
independent invertible matrices. Then, the solution PmiO = YlT=i Xi{i)Pi of (|^ is such that 
the vector A(.^) G M”* satishes M^{^)X{^) = with 

= trace(H^H^(e)i^^ PjA{0V) and (0 = trace(H'^P,H(0H). (6) 

The procedure for the computation of and is given in Algorithm]^ Note 

that only mK matrix-vector products involving the implicit matrices P* are required. 

Now the question is to choose a matrix V such that ||(/ — PA{^))V\\f provides a good 
approximation of ||/ — PA(^)||j7’ for any P eY^ and ^ G S. 
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Algorithm 1 Computation of M^{^) and S^{^) 

Require: A(^), {Pi,, Pm} and V = {vi,... ,vk) 

Ensure: M^(^) and S^{^) 

1: Compute the vectors Wi^k = Pi^iO'^k ^ for 1 < A; < A and 1 < i < m 
2: Set Wi = ..., Wi^x) G 1 < -i < m 

3: Compute = trace(hh7’iyj) for I < i,j < m 

4: Compute Sy= trace(C^lEi) for 1 < i < m 


2.2.1 Hadamard matrices for the estimation of the Probenius norm of an im¬ 
plicit matrix 

Let B an implicit u-by-n matrix (consider B = I—PA{^), with P G Ym and G S). Following 
[3] , we show how Hadamard matrices can be used for the estimation of the Frobenius norm of 
an implicit matrix. The goal is to hnd a matrix V such that ||PH||i7’ is a good approximation 
of ||P||f- The relation ||P1/|||, = tia.ce{B'^BVV'^) suggests that V should be such that VV'^ 
is as close as possible to the identity matrix. For example, we would like V to minimize 


err(Vy 


1 

n(n — 1) 




i=l j^i 


l-vv^\\l 

n{n — 1) ’ 


which is the mean square magnitude of the off-diagonal entries of VV'^. The bound err{V) > 
\/{n — K)/{{n — 1)A) is known to hold for any V G whose rows have unit norm |39]. 

Hadamard matrices can be used to construct matrices V such that the corresponding error 
erriy) is close to the bound, see [3]. 

A Hadamard matrix Hg is a s-by-s matrix whose entries are ±1, and which satisfies 
HgHj = si where I is the identity matrix of size s. For example. 


H,= 



1 

-1 


is a Hadamard matrix of size s = 2. The Kronecker product (denoted by (g)) of two Hadamard 
matrices is again a Hadamard matrix. Then it is possible to build a Hadamard matrix whose 
size s is a power of 2 using a recursive procedure: H 2 k+i = H 2 ® ^ 2 '=- The (i, j)-entry of 
this matrix is (—1)“^^, where a and b are the binary vectors such that i = X]fc>o 
j = For a sufficiently large s = 2^ > max(n. A), we define the rescaled partial 

Hadamard matrix V G as the first n rows and the first A columns of Hg/y/K. 
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2.2.2 Statistical estimation of the Frobenius norm of an implicit matrix 

For the computation of the Frobenius norm of S, we can also use a statistical estimator as 
hrst proposed in [26]. The idea is to dehne a random matrix V G with a suitable 

distribution law V such that ||i?ld||i7’ provides a controlled approximation of ||i?||F with high 
probability. 

Definition 2.2 A distribution V overW^^^ satisfies the {e, 5)-concentration property if for 
all B e 

P(|||BV'||5.-||B||J,|>E||B||p<i. (7) 


Two distributions T> will be considered here. 

(a) The rescaled Rademacher distribution. Here the entries of H G are independent 

and identically distributed with Vij = with probability 1/2. According to Theorem 

13 in |2], the rescaled Rademacher distribution satisfies the (e:, 5)-concentration property for 

K >Qe-Hn{2n/6). (8) 


(b) The subsampled Randomized Hadamard Transform distribution (SRHT), first intro¬ 
duced in [1]. Here we assume that u is a power of 2. It is defined by R = G 

^nxK 


• D G is a diagonal random matrix where Dj j are independent Rademacher random 
variables (i.e. Di^i = ±1 with probability 1/2), 






G 


is a Hadamard matrix of size n (see Section 2.2.1), 


R G is a subset of K rows from the identity matrix of size n chosen uniformly 

at random and without replacement. 


In other words, we randomly select K rows of Hn without replacement, and we multiply the 
columns by We can find in [371|7| an analysis of the SRHT matrix properties. In the 

case where n is not a power of 2, we define the partial SRHT (P-SRHT) matrix V G 
as the first n rows of a SRHT matrix of size s x K, where s = jg the smallest power 

of 2 such that n < s < 2n. The following proposition shows that the (P-SRHT) distribution 
satisfies the {e, 5)-concentration property. 


Proposition 2.3 The (P-SRHT) distribution satisfies the {e, 5)-concentration property for 

K > 2{e^ - ln(4/5)(l + ^8 \n{An/5)f. (9) 
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Proof: Let B G We define the square matrix B of size s = ^ whose hrst 

n X n diagonal block is B, and 0 elsewhere. Then we have ||i?||F = ||-B||f- The rest of the 
proof is similar to the one of Lemma 4.10 in [^. We consider the events A = {(1—e)||i?||p < 
< (1 + e:)||i?|||.} and E = {max* \\BDHj'ei\\l ^ (1 + \/8ln(2s/(5))^||i?|||.}, where e* 
is the Tth canonical vector of The relation P(y4'^) < P(y4'^|ii^) + P(£'‘^) holds. Thanks 
to Lemma 4.6 in [7] (with t = a/ 8 ln(2s/(5)) we have P(£^/ < 5/2. Now, using the scalar 
Chernoff bound (Theorem 2.2 in |37j with fc = 1) we have 

F{A^\E) = n\\BV\\l < (1 - e)\\B\\l\E) + P(|| W||| > (1 + e)\\B\\l\E) 

< (e“^(l — ln(2s/5))-2 (^l + g-^-l-£^-f<'(l+A/81n(2s/5))-2 

< 2(e'^(l + ln(2s/<5))-2^ 

The condition ([^ implies ¥{A^\E) < 5/2, and then P(74'^) < 5/2 + 5/2 = 5, which ends the 
proof. ■ 


Such statistical estimators are particularly interesting for that they provide approxima¬ 
tions of the Frobenius norm of large matrices, with a number of columns K for V which scales 
as the logarithm of n, see (|^ and ([^. However, the concentration property ([^ holds only 
for a given matrix B. The following proposition 2H extends these concentration results for 
any matrix i? in a given subspace. The proof is inspired from the one of Theorem 6 in [15] . 
The essential ingredient is the existence of an £-net for the unit ball of a hnite dimensional 
space (see 0 )- 


Proposition 2.4 Let V G be a random matrix whose distribution T> satisfies the 

(e, 5)-concentration property, with £ < 1. Then, for any L-dimensional subspace of matrices 
Ml C and for any C > 1, we have 

PdllBrllJ, - ||B||J.| > t(C+ 1)/(C - 1)||B||J,,VB 6 Mt) < (9C/e)H. (10) 

Proof: We consider the unit ball Bl = {B G Ml ■ ||i?||F < 1} of the subspace Ml- It is 
shown in |B] that for any s' > 0, there exists a net Aff C Bl of cardinality lower than ifi/'s)^ 
such that 

min_ \\B — B-g\\p < e, \/B G Bl- 

In other words, any element of the unit ball Bl can be approximated by an element of N[ 
with an error less than s'. Using the {e, 5)-concentration property and a union bound, we 
obtain 

\\\BA^\\l - \m\l\ < e\\B4l, WB^eXl 


( 11 ) 









with a probability at least 1 — 5(3/e)'^. We now impose the relation e = e:/(3C), where 
C* > 1. To prove (10), it remains to show that equation ([TT|) implies 


IIIWII 


mil 


<e{C + l)/{C-l)\\B\\l, WBgMl. 


( 12 ) 


We dehne B* G argmax^gg^ |||i?l/|||, —||i?|||,|. Let B^ G Ml be such that \\B* —B^\\p < e, 
and B~ = argminggspan(Bj) \\B* — B\\p. Then we have \\B* — B~\\p = ||-B*|||' — ||-Si|||, < 
and {B* — B~,B~) = 0, where (•, •) is the inner product associated to the Frobenius norm 
II ■ II i?. We have 


r] := |||5W|||- ||S*|||| = \\\{B* - B^)V + B*~Vfp - \\B* - B*~ + B*~\\],\ 

= |||(S* - B^)Vfp + 2{{B* - B*~)V,B^V) + \\B^Vfp - \\B* - B*~\\l - \\B^\\l\ 

< IIKS* - BI)V\\l -\\B*- Bi\\l\ + \mv\\l - \m\\l\ + MB* - B*~)V\\p\\B*~V\\p 

We now have to bound the three terms in the previous expression. Firstly, since {B* 
Bi 


\B* - B^\\p G Bl, the relation |||(5* - Bi)V\\l - \\B* - B^\\j,\ < \\B* - 5i|||r/ < Pr] 

< e. Thirdly, by dehnition of p, we 


\B 


* 112 


< e\\B. 


* 112 


holds. Secondly, (11) gives \ \\B~V\\\ — ^ c||jj~||g 

can write \\{B* - B*~)V\\l < {I + 7])\\B* - B*~\\l <^{1 + r]) and \\B*~Vfp < (l + £)||5i||^ < 
1 +£, so that we obtain 2|| {B* — i?i)K||ir||i?~l/||g < 2e\/lPl:y/Y^Pri. Finally, from (11), we 
obtain 

rj <Pri + e + 2£Vl + £\/l + V (13) 


Since £ < 1, we have £ = e/{?>C) < 1/3. Then P < e and \flP~e < 3/2, so that (13) implies 
Tj <erj + e + 3e\/l + rj <eri + e + 3£(1 + p/2) < Serj + e + 3£, 


and then p < (e + 3£)/(l — 3P) < e(C + 1)/(C — 1). By dehnition of p, we can write 

IIIWII) 


II 


l^llll 


< e{C + 1)/(C — 1) for any B G Bp, that implies (12). 


Proposition 2.5 Let G S, and let PmiO ^ ^ defined by (§ where V G is a 

realization of a rescaled Rademacher matrix with 

K > 6e-^ ln(2n(9C'/e)"^+Vh), (14) 

or a realization of a P-SRHT matrix with 

K > 2{P - e^S)-^ ln(4(9C'/£)™+V(5)(l + ^8 \n{An{fiC/6)^^+^/5)f (15) 

for some 5 > 0, £ < 1 and C > 1. Assuming s' = e{C + 1)/(C' — 1) < 1, 

||/ - ^’„K)^K)IIf < Jpf, min ||/ - FA(0 ||f (16) 

V -L C 

holds with a probability higher than 1 — h. 
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Proof: Let us introduce the subspace M^+i = YmA{^) + span(/) of dimension less than 
m + 1, such that {/ — PA{^) : P G Ym} C Mm+i- Then, we note that with the conditions 
(14) or (15), the distribution law P of the random matrix V satishes the (e, 
concentration property. Thanks to Proposition 2A the probability that 

Ii|(/ - P41(0)P||| - 11/ - PAiOWll < e'\\I - PAiOWl 

holds for any P G Ym is higher than 1 — 5. Then, by dehnition of Pm{0 have with a 

probability at least 1 — 5 that for any P G Y^, it holds 


|/- -PmmiOMF, 


< 


^^||(/ - P/l(f))V'||^ < PS||/ - P.4(f)|lF. 


Then, taking the minimum over P G Y^, we obtain (16). 


Similarly to Proposition |2^ we obtain the following properties for with Pm(0 

the solution of ([^. 


(17) 


(18) 


Proposition 2.6 Under the assumptions of Proposition [^75| the inequalities 

(1 - QV„«))"(1 - f')-' < II(^ - ft.K).4(f))r|i5. < n (1 - (1 - F')a^(f)) 

and 

«(Pm(0-4(0) < - (u - l)a^(0 

hold with probability 1 — 5, where am{i) is the lowest singular value of Pm{i)A{ff). 

Proof: The optimality condition for Pm(0 yields ||(J—Pm(0^(0)^ill’ = ll^ll' 

Since Pm(0^(0 ^ ^m+i (where M^+i is the subspace introduced in the proof of Proposition 
(2.5)), we have 

WPmmmwi > (1 - e')\\Pmmmi ( 19 ) 

with a probability higher than 1 — 5. Using i|U|||’ = n (which is satishes for any realization of 
the rescaled Rademacher or the P-SRHT distribution), we obtain ||(/ — Pm{^)A{^))V\\‘jp < 
n — (1 — P)||Pm(0^(0llF with a probability higher than 1 — 5. Then, ||Pm(0^(0111’ — 
naraiO'^ yields the right inequality of ( pT| . Following the proof of Lemma 2.6 in [21], we 
have (1 — Om(O^) — II/ ~ /m(0^(0 IIf- Together with (19), it yields the left inequality of 
Furthermore, with probability 1 — 5, we have n — (1 — P)||Pm(0^(0llF — 0- Since 


iPmiomvwi. 
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the square of the Frobenius norm of matrix is the sum of squares of its singular 

values, we deduce 


(n - l)aUO^ + < ||P^(0A(0lll < ^(1 - 

with a probability higher than 1 — 5, where /3m{0 i® i'h® largest singular value of Pm(^)A(^). 
Then (18) follows from the dehnition of K{Pm{0^i0) = f^miO /■ 


2.2.3 Comparison and comments 


We have presented different possibilities for the dehnition of V. The rescaled partial Hadamard 


matrices introduced in section 2.2.1 have the advantage that the error err{V) is close to the 


theoretical bound \/{n — K)/{{n — 1)P), see Figure 1(a) (note that the rows of V have 
unit norm). Furthermore, an interesting property is that VV^ has a structured pattern 


(see Figure 1(b)). As noticed in [3], when K = 2^ the matrix VV'^ have non-zero entries 
only on the 2'^^-th upper and lower diagonals, with A; > 0. As a consequence, the error on 
the estimation of ||P||f will be induced only by the non-zero off-diagonal entries of B that 
occupy the 2'^^-th upper and lower diagonals, with fc > 1. If the entries of B vanish away 
from the diagonal, the Frobenius norm is expected to be accurately estimated. Note that 
the P-SRHT matrices can be interpreted as a “randomized version” of the rescaled partial 


Hadamard matrices, and Figure 1(a) shows that the error err{V) associated to the P-SRHT 
matrix behaves almost like the rescaled partial Hadamard matrix. Also, P-SRHT matrices 
yield a structured pattern for VV'^, see Figure 1(c) The rescaled Rademacher matrices give 


higher errors err{V) and yield matrices VV^ with no specihc patterns, see Figure 1(d) 


The advantage of using rescaled Rademacher matrices or P-SRHT matrices is that the 
quality of the resulting projection Pm(0 can be controlled with high probability, pro¬ 
vided that V has a sufficiently large number of rows K (see Proposition 2.5). Table 


shows the theoretical value for K in order to obtain the quasi-optimality result (16) with 


+ £')/(! ~ ^0 = 10 ^^0 5 = 0.1%. It can be observed that K grows very slowly with 
the matrix size n. Also, K depends on m linearly for the rescaled Rademacher matrices and 


quadratically for the P-SRHT matrices (see equations (14) and (15)). However, these theo¬ 
retical bounds for K are very pessimistic, especially for the P-SRHT matrices. In practice, 
it can be observed that a very small value for K may provide very good results (see Section 
1^. Also, it is worth mentioning that our numerical experiments do not reveal signihcant dif¬ 
ferences between the rescaled partial Hadamard, the rescaled Rademacher and the P-SRHT 
matrices. 
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(c) Distribution of the entries of (in 

absolute value) where y is a sample of the 
P-SRHT matrix with K = 100. 



(b) Distribution of the entries of VV'^ (in 
absolute value) where V is the rescaled par¬ 
tial Hadamard matrix with K = 100. 


too 


200 


300 


400 


500 


600 


0.2 to 0.25 (8344 entries) 
3.25 to 0.3 (3434 entries) 
1.3 to 0.99 (626 entries) 

■ 0.99 to 1 (600 entries) 



too 200 300 400 


500 600 


(d) Distribution of the entries of VV'^ (in 
absolute value) where P is a sample of the 
rescaled Rademacher matrix with K = 100. 


Figure 1: Comparison between the rescaled partial Hadamard, the rescaled Rademacher and 
the P-SRHT matrix for the dehnition of matrix V, with n = 600. 


2.3 Ensuring the invertibility of the preconditioner for positive 
definite matrix 

Here, we propose a modihcation of the interpolation which ensures that Pm{0 is invertible 
when H(,^) is positive dehnite. 

Since is positive dehnite. Pi = is positive dehnite. We introduce the vectors 
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(a) Rescaled Rademacher distribution. 



m = 2 

m = 5 

m = 10 

m = 20 

TO = 50 

n = 10"^ 

239 

363 

567 

972 

2 185 

n = 10® 

270 

395 

599 

1 005 

2 219 

n = 10® 

301 

427 

632 

1 038 

2 253 


(b) P-SRHT distribution. 



TO = 2 

TO = 5 

TO = 10 

TO = 20 

TO = 50 

n = 10^ 

27 059 

63 298 

155 129 

455 851 

2 286 645 

n = 10® 

30 597 

69 129 

164 750 

473 011 

2 326 301 

n = 10® 

34 112 

74 929 

174 333 

490 126 

2 365 914 


Table 1: Theoretical number of columns K for the random matrix V in order to ensure (16), 
with ■\/(l + £')/(I — £') = 10 and 6 = 0.1%. The constant C has been chosen in order to 
minimize K. 


7“ e 


and 7 ^ G whose components 


. {PiW,w) ^ j , {PiW,w) 

7,- = mt — r —^7-— > 0 and 7, = sup —n—— < 00 

correspond respectively to the lowest and highest eigenvalues of the symmetric part of Pj. 
Then, for any P = 1 AjP^ G Y^, 


, {Pw,w) ,. 

\i 112 ^7) - (-^ ,7 ), 


( 20 ) 


where > 0 and A“ > 0 are respectively the positive and negative parts of A = A"*" — 
A“ G M™'. As a consequence, if the right hand side of (20) is strictly positive, then P is 
invertible. Furthermore, we have ||P|| < (A’*' + X~,C), where C G M”* is the vector of 
component Cj = ||Pj||, where ||Pj|| denotes the operator norm of Pj. If we assume that 
(A+, 7 “) — (A“, 7 +) > 0, the condition number of P satisfies 


k(P) = ||P| 


>-ii 


< IIPI 


inf 


{Pw, w) 


\w\ 


-1 


< 


(A+ + A-,C') 


(A+, 7 -) - (A-, 7 +)' 

It is then possible to bound k(P) by R by imposing 

(A+ + A-, C) < /t((A+, 7 -) - (A-, 7 +)), 

which is a linear inequality constraint on A"*" and A“. We introduce two convex subsets of 
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Ym defined by 


A+ > 0, A- > 0 

= -! -Y^iP ■ p 0 

(A+, Rj- -C)- (A-,K7+ + C)>0 


i=l 


2 = 1 


= E pp ■ pp'^ 


2=1 


From (20), we have that any nonzero element of is invertible, while any nonzero 
element of Y^ is invertible and has a condition number lower than R. Under the condition 
R > maXiCj/yj”, we have 

Y* C Y‘ c (22) 

Then dehnitions ([^ and (|^ for the approximation Pm(0 can be replaced respectively by 


PmiO = argmin ||/- P24(^)||i7’, (23a) 

PeY+ or 

PmiO = argmin \\{I - PA{^))V\\f, (23b) 

P£Y+ or Y^ 

which are quadratic optimization problems with linear inequality constraints. Furthermore, 
since Pi G Y^ for all i, all the resulting projections Pm(0 interpolate 24(^)“^ at the points 

^li ■ ■ ■ 1 ^m- 

The following proposition shows that properties ([^ and (|^ still hold for the precondi¬ 
tioned operator. 


Proposition 2.7 The solution Pm{C) of (23a) is such that satisfies (|^ and 0. 


Also, under the assumptions of Proposition 2.5, the solution Pm{0 of (23b) is such that 
PmiO^iO satisfies im and ([I^ with a probability higher than 1 — 5. 


Proof: Since Yfi^ (or Yfif) is a closed and convex positive cone, the solution PmiO of (23a) 
is such that trace{{I — PmiO^iOViPmiO — P 0 fo^ P ^ (or Ym)- Taking 

P = 2Pm{0 ^-ricl P = 0, we obtain that trace((/ — Pm{0^iO)^PmiOP^iO) = 0; which im¬ 
plies ||Pm(0^(0llF = trace(Pm(^)y4(^)). We refer to the proof of Lemma 2.6 and Theorem 
3.2 in to deduce ([^ and (|^. Using the same arguments, we prove that the solution 
PmiO of (23b) satishes ||Pm(0^(0f^llF = trace(U^Pm(0^(0f^)) ^rid then that ( p(^ and 


(18) hold with a probability higher than 1 — 5. 
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2.4 Practical computation of the projection 

Here, we detail how to efficiently compute M^{^) and S^{^) given in equation ([^ in a multi¬ 
query context, i.e. for several different values of The same methodology can be applied 
for computing M(^) and S{^). We assume that the operator has an affine expansion 
of the form 

niA 

= (24) 

k=l 


where the Ak are matrices in and the : 

M^{^) and also have the affine expansions 


-)■ 


are real-valued functions. Then 


ruA rriA 

"m(«) = E E tv^ceiV-^AlPlP,A,V), 

k=l 1=1 


SYiO = > >fc(0 trace(W'P,HfcH), 


mA 


k=l 


(25a) 


(25b) 


respectively. Computing the multiple terms of these expansions would require many com¬ 
putations of traces of implicit matrices and also, it would require the computation of the 
affine expansion of 4(^). Here, we use the methodology introduced in [9] for obtaining affine 
decompositions with a lower number of terms. These decompositions only require the knowl¬ 


edge of functions in the affine decomposition (24), and evaluations of Mh(^) and SY (0 
(that means evaluations of 24(,^)) at some selected points. We briefly recall this methodology. 

Suppose that g : E ^ X, with X a vector space, has an affine decomposition g{^) = 
YYk=iCk{0gk, with (Cfc : S —)■ M and gu G X. We first compute an interpolation of C(^) = 
(Ci(0, ■ ■ ■ ,Cm(0) under the form ((0 = ^fc(0C(O, with mg < m, where 

are interpolation points and \['i(^),..., (0 fh® associated interpolation functions. Such 


an interpolation can be computed with the Empirical Interpolation Method [29] described 
in Algorithm]^ Then, we obtain an affine decomposition g(^) = YYh ^kiCldiU) which can 
be computed from evaluations of g at interpolation points Yl- 

Applying the above procedure to both and we obtain 


PPM 


ms 


M''(o « 4-^(0 MpjB, s''(o « s'-ya. 


(26) 


k=l 


k=l 


The hrst (so-called offline) step consists in computing the interpolation functions "^kiO 
and 4/fc(^) and associated interpolation points and Q using Algorithm with input 
{^i^j}i<i,j<mA und {^i}i<i<mA respectively, and then in computing matrices {Yk) und 
vectors {ff) using Algorithm hi The second (so-called online) step simply consists in 


computing the matrix M^{^) and the vector S'^(^) for a given value of ^ using (26). 
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Algorithm 2 Empirical Interpolation Method (EIM). 

Require: (Ci(-), • • • , Cm(-)) 

Ensure: d'i(-), ..., d'fc(-) and ..., 

1: Define = (i{^) for all 

2: Initialize e = 1, fc = 0 

3: while e > tolerance (in practice the machine precision) do 
4: k = k + 1 

5: Find {il,Ck) ^ argmax |Rfc(*,OI 

6: Set the error to e = \Rk{il^^k)\ 

7: Actualize Rk+i{i,0 = Rkii,0 - Rk{i,il)Rk{il,i)/Rk{il,il) for all 

8: end while 

9: Fill in the k-hy-k matrix Q : Qjj = for all ^ < i-, j < k 

10: Compute for all ^ and 1 <i <k 


3 Preconditioners for projection-based model reduc¬ 
tion 


We consider a parameter-dependent linear equation 

/!«)««) = 6 «). 


(27) 


with A(^) G and h{^) = M”. Projection-based model reduction consists in projecting 

the solution u{^) onto a well chosen approximation space C X := M” of low dimension 


r <^n. Such projections are usually defined by imposing the residual of (27) to be orthogonal 


to a so-called test space of dimension r. The quality of the projection on X^ depends on the 
choice of the test space. The latter can be defined as the approximation space itself X^, thus 
yielding the classical Galerkin projection. However when the operator A(^) is ill-conditioned 
(for example when A(,^) corresponds to the discretization of non coercive or weakly coercive 
operators), this choice may lead to projections that are far from optimal. Choosing the test 
space as {R'^A{^)vr : G X^}, where R'^A{^) is called the “supremizer operator” (see 

e.g. [35]), corresponds to a minimal residual approach, which may also results in projections 
that are far from optimal. In this section, we show how the preconditioner Pm{0 can be 
used for the definition of the test space. We also show how it can improve the quality of 
residual-based error estimates, which is a key ingredient for the construction of suitable 
approximation space X^ in the context of the Reduced Basis method. 

X is endowed with the norm || ■ ||x defined by || ■ ||^ = {Rx-, ■), where Rx is a symmetric 
positive definite matrix and {•, •) is the canonical inner product of M"". We also introduce the 
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dual norm || ■ ||x' = ' il^ such that for any v,w G X we have |(u,ta)| < ||u||x Utr'Hx'- 


3.1 Projection of the solution on a given reduced subspace 

Here, we suppose that the approximation space Xr has been computed by some model order 
reduction method. The best approximation of u{^) on Xr is u*{^) = argmin^igx,. \\u{^) — "ylU 
and is characterized by the orthogonality condition 

RxVr) = 0, WVr E Xr, (28) 


or equivalently by the Petrov-Galerkin orthogonality condition 


{A(0<(?) - 6(^), A-^{0Rxv,} = 0, Vt., e X,. 


(29) 


Obviously the computation of test functions A~'^{^)RxVr for basis functions Vr of Xr is 
prohibitive. By replacing by PmiO^ obtain the feasible Petrov-Galerkin formula¬ 

tion 


- HO: P^iORxVr) = 0 , yvr E Xr- 


Denoting by G G a matrix whose range is Xr, the solution of (30) is Ur{0 
where the vector a{0 E is the solution of 


(30) 

UaiO 


{U^RxPmiOAm)a{0 = {U^RxPmiOHO)- 


Note that (30) corresponds to the standard Galerkin projection when replacing PmiO 


by Rj)^ . Indeed, the orthogonality condition (30) becomes (24(^)Mr(0 ~ HOXr) = 0 for all 
Vr E Xr- 


Remark 3.1 Here, the preconditioner Pm{0 for the definition of the parameter- 

dependent test space {PfiiififiRxVr '■ Vr E Xr} which defines the Petrov-Galerkin projection 
(30). However, Pmi.0 could also be used to construct preconditioners for the solution of the 


linear system (U"^A{0U)a{0 = {U'^b{0) corresponding to the Galerkin projection on Xr- 
Following the idea proposed in m, such preconditoner can take the form (U'^Pm{^)U), thus 
yielding the preconditioned reduced linear system 

{Wp„,(0U){WA(()u)a(0 = {u^p„(t;)u){u'^HO)- 

Such preconditioning strategy can be used to accelerate the solution of the reduced system of 


eguations when using iterative methods. However, and contrarily to (30), this strategy does 


not change the definition ofUr{^), which is the standard Galerkin projection. 
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We give now a quasi-optimality result for the approximation This analysis relies 

on the notion of h-proximality introduced in [T3] . 

Proposition 3.2 Let Sr,m(0 ^ [O 5 1] defined by 

\Vr - Rx\Pm{OMOfRxWr\\x 


dr,m(0 — max min 


\Vr\\x 


The solutions u*(^) E and Mr(0 ^ of (28) and (30) satisfy 

|K(0 - w(0llx < - w(0IU- 

Moreover, if < 1 holds, then 

ll«K) - «,K)IU < (1 - (?) - <(?)IU. 


Proof: The orthogonality condition (28) yields 


«(0 -UriO^RxVr) = ( m (0 -Ur{f),RxVr) = (&(0 '^{ORxVr) 


for all Vr E Xr- Using (30), we have that for any Wr E W, 

«(0 - Ur{0,RxVr) = m) - A{0Ur{0,A-^i0RxVr " Pm{0^ RxWr), 

= {u{0 - UriO: RxVr - {Pm{0M0)'^RxWr), 

< l|w(e) - W.(0IU WRxVr - {PmiOMOfRxWrWx’ 
= I|W(0 - Wr-(0IU ll^^r- - Rx\Pm{0M0fRxWr\\x. 
Taking the inhmum over Wr E Xr and by the dehnition of 6r,m{0j obtain 

{u*{f) - Ur{0,RxVr) < 5r,m(0 ll«(0 “ W(0 lU ll^'rUx- 
Then, noting that u*(f) — Ur(f) E Xr, we obtain 

II */'C) /'Oil — {'0‘ri^) ~ 'Orif): Rx^r) , ^ /^nii /■ r\ 

Vr&Xr |Tr||X 


that is (32). Finally, using orthogonality condition ( |28[ ), we have that 

ll«(?) - «r(?)||i = ||«(?) - <(?)||i + ll<(?) - «r(?)||i, 

< ii«(?) - «:(?)iit+v™(?)^ii“(?) - “r(?)iiy 


(31) 

(32) 

(33) 


from which we deduce (33) when hr,m(0 < 1- 


An immediate consequence of Proposition 3.2 is that when 5r,m(0 = 0; Petrov- 
Galerkin projection Mr(0 coincides with the orthogonal projection Following [H], we 

show in the following proposition that (5r,m(0 can be computed by solving an eigenvalue 
problem of size r. 
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Proposition 3.3 We have Sr,m(0 = \/l — 7 , where 7 is the lowest eigenvalue of the gener¬ 
alized eigenvalue problem Cx = '^Dx, with 

C = e 

D = U^RxU e 


where B = {Pm{^)A{f,))'^RxU G and where U G is a matrix whose range is Xr- 
Proof: Since the range of U is we have 

aeK’-beR’" ||f^a||x 

For any a E W, the minimizer b* of ||?7a — Rx^BbWj^ over 6 G M'" is given by b* = 
{B^RflB)~^B^Ua. Therefore, we have \\Ua — RflBb*\W = ||f/a|||^ — {Ua,Bb*), and 

{U^B{B^R-^B)-^B^Ua, a) 


= 1 - inf 

aSR'- 


{U^RxUa, a) 


which concludes the proof. 


3.2 Greedy construction of the solution reduced subspace 

Following the idea of the Reduced Basis method [36l|3H], a sequence of nested approximation 
spaces {Xr}r>i in X can be constructed by a greedy algorithm such that X^+i = Xr + 
span(-u(,^,?,^)), where is a point where the error of approximation of u{f) in Xr is 

maximal. An ideal greedy algorithm using the best approximation in Xr and an exact 
evaluation of the projection error is such that 


uKf) is the orthogonal projection of u{f) on Xr defined by (28), 


fr+i e argmax \\u{f) -<(Ollx- 


CeH 


(34a) 

(34b) 


This ideal greedy algorithm is not feasible in practice since u{f) is not known. Therefore, 
we rather rely on a feasible weak greedy algorithm such that 


Ur{^) is the Petrov-Galerkin projection of u{f) on Xr defined by (30 


cRB 

•sr+l 


argmax ||Pm( 0 (A( 0 «r -(0 “ ^( 0 )IU. 


«6S 


Assume that 


(35a) 

(35b) 


\\u{0 - «r-(0IU < ll^m(0(^(0«r-(0 “ K0)\\x < f^m ||m(0 “ «r-(0IU 
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holds with = inf^g=am(0 > 0 and /3m = supgg=/3m(0 < o®; where 0^(0 PmiO 
respectively the lowest and largest singular values of Pm (0^(0 with respect to the norm 
II ■ ||x, respectively dehned by the inhmum and supremum of ||Pm(0^(0'^ll^ v G X 


such that ||n||x = 1- Then, we easily prove that algorithm (35) is such that 


l“Kr+l) -“rK,+l)|Lv > 7mm|gt||uK) -MrK)||x, 


(36) 


where 7m = a^/ /3m < 1 measures how far the selection of the new point is from the ideal 


greedy selection. Under condition (36), convergence results for this weak greedy algorithm 
can be found in 0 US] ■ 

We give now sharper bounds for the preconditoned residual norm that exploits the fact 
that the approximation Ur{^) is the Petrov-Galerkin projection. 


Proposition 3.4 LetUr{i) the Petrov-Galerkin projection ofu{^) on defined by (29) 
Then we have 


ll«(0 -«r(OIU < -KO)IU < ll«(0 - w(OIU, 


with 


ar,m{0 = inf sup 

’'SX Wr&X, 


I3r,mi0 = snp inf 

vGX '^r^Xr 


wiPUQAiorRxvWx' ^ 

\\V-Wr\\x 

\\{Pmm{OrRx{v-Wr)\\x' 

l|n|U 


Proof: For any v G X and Wr € Xr and according to (30), we have 

(w(0 - UriO, Rxv) = m) - A-^iORxV - PmiORxWr) 

= (PmiOm) - AiOUriO), {PmiOMOr^RxV - RxWr) 

< \\R\\x \\{Pm{0M0)~^RxV - RxWrWx', 

where P(0 •= P-m/OiKO ~ ■ Taking the inhmum over Wr G Xr, dividing by ||n||x 

and taking the supremum over n G X, we obtain 


|m( 0 -Mr(0llx < ||P(0llxsup inf 

vGX 'R’r^Xr 

= ll^(0IUsup inf 


||(Pm(0^(e))-^^x^^-^xi^.||x' 


ITIU 

\V - Wr\\x 


= ll^(0iU ( inf sup 

^veX rOr&Xr 


„^^X->r^Xr\\{Pm{OmYRxv\\x'' 

\\{PmmY)f RXVWX' 


|n - Wr\\x 


-1 
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which proves the first inequality. Furthermore, for any u G X and Wr G X^, we have 


{P^{Om)-mUr{i)),RxV) = m-mnr{O.Pl{i)Rx{v-Wr)) 

< ll“(0 “ \\{Pm{OMO)'^Rx{v - Wr)\\x'- 

Taking the infimum over Wr G X^, dividing by ||n||x and taking the supremum over n G X, 
we obtain the second inequality. ■ 


Since d X^_|_]^, we have ttr+l,m(0 — ^ /^r,m(0 — /^m(0' 

Equation (36) holds with 7 ^ replaced by the parameter 7 ^,^ = a^m/ Pr,m- Since 7 ,.,^ 
increases with r, a reasonable expectation is that the convergence properties of the weak 
greedy algorithm will improve when r increases. 


Remark 3.5 When replacing PmiO Rx^j ^^e preconditioned residual norm ||Pm(0(^(0'^r-(0“ 
^(0)llx turns out to be the residual norm ||74(,^)mj,(0 — which is a standard choice 

in the Reduced Basis method for the greedy selection of points (with Rx being associated 
with the natural norm on X or with a norm associated with the operator at some nom¬ 
inal parameter value). This can be interpreted as a basic preconditioning method with a 
parameter-independent preconditioner. 


4 Selection of the interpolation points 

In this section, we propose strategies for the adaptive selection of the interpolation points. 
For a given set of interpolation points ^ 1 ,... ,^m, three different methods are proposed for 
the selection of a new interpolation point The first method aims at reducing uniformly 

the error between the inverse operator and its interpolation. The resulting interpolation of 
the inverse is pertinent for preconditioning iterative solvers or estimating errors based on 
preconditioned residuals. The second method aims at improving Petrov-Galerkin projections 
of the solution of a parameter-dependent equation on a given approximation space. The 
third method aims at reducing the cost for the computation of the preconditioner by reusing 
operators computed when solving samples of a parameter-dependent equation. 

4.1 Greedy approximation of the inverse of a parameter-dependent 
matrix 

A natural idea is to select a new interpolation point where the preconditioner Pm{0 is not 
a good approximation of y4(,^)“^. Obviously, an ideal strategy for preconditioning would be 
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to choose ^m+i where the condition number of is maximal. The computation of 

the condition number for many values of ^ being computationaly expensive, one could use 
upper bounds of this condition number, e.g. computed using SCM izg. 

Here, we propose the following selection method: given an approximation Pm{0 associ¬ 
ated with interpolation points ^i,..., a new point .^m+i is selected such that 


e argmax ||(J - 

«6S 


(38) 


where the matrix V is either the random rescaled Rademacher matrix, or the P-SRHT matrix 
(see Section 2.2). This adaptive selection of the interpolation points yields the construction 
of an increasingsequence of subspaces Ym+i = + span(H(^m+i)~^) in R = This 

algorithm is detailed below. 


Algorithm 3 Greedy selection of interpolation points. 


Require: M. 

Ensure: Interpolation points ^i,... and interpolation Pm{0- 
1 : Initialize Po(0 = I 
2 : for m = 0 to M — 1 do 

3: Compute the new point according to (38) 

Compute a factorization of A(^m+i) 

Dehne as an implicit operator 


4: 

5: 

6 : 

7: 


Update the space Rm+i = Ym + span(A(,^m-i-i) 
Compute Pm+iiO = argminpgy^^^ \\{I - PA{^))V\\f 

end for 


The following lemma interprets the above construction as a weak greedy algorithm. 


Lemma 4.1 Assume that A(,^) satisfies agll ' II ^ 11^(0 ‘ II < /^oll ■ || for all ^ G S, and let 
PmiO be defined by (|^. Under the assumption that there exists e E [0,1[ such that 

Ilia - pAmvwi -\\i- pmwii <41- pmwi m 


holds for all ^ eE and P E Y^, we have 


\-ii 


\\Pm{fm+i) - Up > y^max min ||P - A(0’ ||f, 

PeYm 


(40) 


with 7 e = Oqa/I — ef [fioy/l + e), and with defined by (38) 


Proof: Since URCUp < ||i?||p||C|| holds for any matrices B and C, with He'll the operator 
norm of C, we have for all P G R, 

iiA(o-' - pia < 11/ - PA(e)iaiiA(o-i < ao4i - pmowf, 

III - PAiOh < ||44(0-^ - PMAm < ^oii/i(o-^ - P\\f- 
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Then, thanks to (39) we have 


PK)-' - p|If < (a,^r^)“'ll(/- -P^rojV'IlF 

and II (/ - Pf1K))1/||f < AVTTillFlK)-' - F||n. 


which implies 


(=||(/ - PFl(f))V'||F < PK)-‘ - PIIf < I - ||(J - Pdl(f))V'||F. 

Po^/l + e 


We easily deduce that .^m+i is such that (40) holds. 


Remark 4.2 The assumption (39) of Lemma f.l can he proved to hold with high probability 
in two cases. A first case is when E is a training set of finite cardinality, where the results of 
Proposition 2^ can be extended to any ^ E E by using a union bound. We then obtain that 
(39) holds with a probability higher than 1 — 5{ffE). A second case is when y4(,^) admits an 
affine decomposition (24) with mA terms. Then the space Ml = span{I — PA{ff) : ^ E E,P E 
Ym} is of dimension L < 1 + mAm and Proposition 2.4 allows to prove that assumption (39) 
holds with high probability. 


The quality of the resulting spaces Wi have to be compared with the Kolmogorov m-width 
of the set y4“^(S) := ^ G S} C K, dehned by 

d,n{A-^{E))Y = min sup mm ||kl(0“^ - P||f, (41) 

c r ? 6 H 

dim(Kn) = m 


which evaluates how well the elements of A ^(S) can be approximated on a m-dimensional 
subspace of matrices. (40) implies that the following results holds (see Corollary 3.3 in [T8]b 


11^(0-' 


^m(OI|F = 


0{m “) 
C>(e-"™') 


if dm(A ^(S))y = 0{m “) 
if dr.{A-\E))y = 0{e-^^'’)' 


where c > 0 is a constant which depends on c and b. That means that if the Kolmogorov 
m-width has an algebraic or exponential convergence rate, then the weak greedy algorithm 
yields an error ||Pm(0 — which has the same type of convergence. Therefore, the 

proposed interpolation method will present good convergence properties when dm{A~^(E))Y 
rapidly decreases with m. 


Remark 4.3 When the parameter set S is [—1,1]'^ (or a product of compact intervals), an 
exponential decay can be obtained when 74(,^)“^ admits an holomorphic extension to a domain 
in C'^ containing S (see m). 
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Remark 4.4 Note that here, there is no constraint on the minimization problem over 
(either optimal subspaces or subspaces constructed by the greedy procedure), so that we have 


no guaranty that the resulting approximations Ym are invertible (see Section 2.3) 


4.2 Selection of points for improving the projection on a reduced 
space 

We here suppose that we want to hnd an approximation of the solution m(^) of a parameter- 


dependent equation (27) onto a low-dimensional approximation space Xr, using a Petrov- 


Galerkin orthogonality condition given by (30). The best approximation is considered as 


the orthogonal projection dehned by (28). The quantity (5r,m(0 dehned by (31) controls 


the quality of the Petrov-Galerkin projection on Xr (see Proposition 3.2). As indicated 
in Proposition 3.3, Sr,m(0 can be efficiently computed. Thus, we propose the following 


selection strategy which aims at improving the quality of the Petrov-Galerkin projection: 
given a preconditioner Pm(0 associated with interpolation points the next point 

^rn+i is selected such that 

e argmax dr,m(0- (42) 

C6S 

The resulting construction is described by Algorithm with the above selection of 
Note that this strategy is closely related with [T3], where the authors propose a greedy 
construction of a parameter-independent test space for Petrov-Galerkin projection, with a 
selection of basis functions based on an error indicator similar to 5r,m(0- 


4.3 Re-use of factorizations of operator’s evaluations - Application 
to reduced basis method 


When using a sample-based approach for solving a parameter-dependent equation (27), the 


linear system is solved for many values of the parameter When using a direct solver for 
solving a linear system for a given a factorization of the operator is usually available and 
can be used for improving a preconditioner for the solution of subsequent linear systems. 

We here describe this idea in the particular context of greedy algorithms for Reduced 
Basis method, where the interpolation points ... ,^r for the interpolation of the inverse 
^(0~^ are taken as the evaluation points ..., for the solution. At iteration r, having 
a preconditioner Pr(0 an approximation Ur(^), a new interpolation point is dehned such 
that 

-RB 




r-\-l 


argmax ||R.(0(A(0Mr(0 “ ^(0 )llx- 

«6H 


Algorithm describes this strategy. 
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Algorithm 4 Reduced Basis method with re-use of operator’s factorizations. 

Require: V, and R. 

Ensure: Approximation Ur{^). 

1 : Initialize Mo(0 = 0; ^ 

2 : for r = 0 to R — 1 do 

3: Find e argmax^gH ||P^(0(A(0Mr(0 “ K0)\\x 

4: Compute a factorization of 

5: Solve the linear system Vr+i = 

6 : Update the approximation subspace A^+i = Xr + span(nr+i) 

7: Define the implicit operator Pr+i = 

8 : Update the space Yr+i (or 

9: Compute the preconditioner : R.+i(^) = argmin ||(/ — PA{^))V\\f 

P&Y,+i{OY 

10 : Compute the Petrov-Galerkin approximation Mr+i(0 of -^r+i using equation 

(0 

11: end for 


5 Numerical results 


5.1 Illustration on a one parameter-dependent model 

In this section we compare the different interpolation methods on the following one parameter- 
dependent advection-diffusion-reaction equation: 


— An -I- v{^) ■ Vm + u = f 


(43) 


defined over a square domain D = [0,1]^ with periodic boundary conditions. The advection 
vector field u(^) is spatially constant and depends on the parameter ^ that takes values in 
[0,1]: u(^) = Dcos(27r,^)ei -|-Dsin(27r^)e2, with D = 50 and ( 61 , 62 ) the canonical basis 
of S denotes a uniform grid of 250 points on [0,1]. The source term / is represented 


in Figure 2(a) We introduce a hnite element approximation space of dimension n = 1600 
with piecewise linear approximations on a regular mesh of D. The mesh Peclet number 
takes moderate values (lower than one), so that a standard Galerkin projection without 
stabilization is here sufficient. The Galerkin projection yields the linear system of equations 
^( 0 w (0 = b, with 


A(,^) = Ao + cos(27r,^)Ai sin(27r,^)A2, 
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where the matrices Ai, A 2 and the vector b are given by 


“1“ 4^i^j 1 (^Ai'ji j / 0j(ei ■ 

Ju Ju 

{A2)i,j = / 0*(e2 ■ V0j) , {b)i = / (j)if, 

Jo. J ^ 


where is the basis of the hnite element 

samples of the solntion. 


space. Figures 2(b), 2(c) and 2(d) show three 



(a) / (b) u(0.05) (c) u(0.2) (d) u(0.8) 

Figure 2: Plot of the source term / (a) and 3 samples of the solution corresponding to 
parameter values ^ = 0.05 (b), .^ = 0.2 (c) and ^ = 0.8 (d) respectively. 


5.1.1 Comparison of the interpolation strategies 


We first choose arbitrarily 3 interpolation points (^1 = 0.05, ^2 = 0.2 and ^3 = 0.8) and show 
the benehts of using the Frobenius norm projection for the dehnition of the preconditioner. 
For the comparison, we consider the Shepard and the nearest neighbor interpolation strate¬ 
gies. Let II ■ ||h denote a norm on the parameter set S. The Shepard interpolation method 
is an inverse weighted distance interpolation: 


HO 


lie- 611 ^^ 

Er=iiie-e, 

1 


if e 7^ 6 
ife = ei 


where s > 0 is a parameter. Here we take s = 2. The nearest neighbor interpolation method 
consists in choosing the value taken by the nearest interpolation point, that means Aj(.^) = 1 
for some i G argmin^ ||^ — ejlls, and Xj{0) = 0 for all j 7 ^ i. 

Concerning the Frobenius norm projection on Ym (or we first construct the affine 
decomposition of M{^) and S{^) as explained in Section 2.4 The interpolation points 


(resp. ^l) given by the EIM procedure for M(^) (resp. S'(^)) are {0.0; 0.25; 0.37; 0.56; 0.80} 
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(resp. {0.0; 0.25; 0.62}). The number of terms niM = 5 in the resulting affine decomposition 
of M(^) (see equation (26)) is less than the expected number m\ = 9 (see equation (25a)). 


Considering the functions $i(0 = *^2(0 = cos(27r,^), <h 3 (.^) = sin(27r.^), and thanks to 

relation cos^ = 1 — sin^, the space 

spanjj{<hj$j} = spanjl, cos, sin, cos sin, cos^, sin^} = spanjl, cos, sin, cos sin, cos^} 
is of dimension ttim = 5. The EIM procedure automatically detects the redundancy in the 


set of functions and reduces the number of terms in the decomposition (26). Then, since the 


dimension n of the discretization space is reasonable, we compute the matrices M{^1) and 


the vectors S{^1) using equation (26). 

The functions Aj(^) are plotted on Figure for the proposed interpolation strategies. 
It is important to note that contrary to the Shepard or the nearest neighbor method, the 
Frobenius norm projection (on or Y^) leads to periodic interpolation functions, i.e. 
Aj(^ = 0) = Ai(^ = 1). This is consistent with the fact that the application ^ i—)■ y4(^) is 
1-periodic. The Frobenius norm projection automatically detects such a feature. 



5 « 


(a) Nearest neighbor 


(b) Shepard (with s = 2) 



Figure 3: Interpolation functions Aj(^) for different interpolation methods. 


Figure 1^ shows the condition number /tm(0 of with respect to We hrst 

note that for the constant preconditioner Pi{^) = ^(^ 2 )”^, the resulting condition number 
is higher than the one of the non preconditioned matrix A(,^) for ^ G [0.55; 0.95]. We 
also note that the interpolation strategies based on the Frobenius norm projection lead to 
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better preconditioners than the Shepard and nearest neighbor interpolation strategies. When 
considering the projection on and (with R = 5x 10"^ such that (22) holds), the resulting 
condition number is roughly the same, so as the interpolation functions of Figures 3(d) and 


3(e) Since the projection on YJ^ requires the expensive computation of the constants 7 +, 7 " 


and C (see Section 2.3), we prefer to simply use the projection on Y^ in order to ensure the 
preconditioner to be invertible. Finally, for this example, it is not necessary to impose any 
constraint since the projection on Y^ leads to the best preconditioner and this preconditioner 
appears to be invertible for any G S. 



e 


-- - Pi(e) = ^-H?2) 

- Nearest neighbor 

Shepard 

- Projection on 17 

- Projection on ITf 

- Projection on Yjjj 


Figure 4: Condition number of Pm(0^(0 for different interpolation strategies. The condi¬ 
tion number of 24(^) is given as a reference. 


5.1.2 Using the Frobenius semi-norm 


We analyze now the interpolation method dehned by the Frobenius semi-norm projection 
Ym ([^ for the different dehnitions of U G proposed in sections 


on 


2 . 2.2 


and 


2 . 2.1 


According to Table the error on the interpolation functions decreases slowly with K 
(roughly as (9(iF“^/^)), and the use of the P-SRHT matrix leads to a slightly lower error. 


The interpolation functions are plotted on Figure 5(a) in the case where K = 8. Even if we 


have an error of 36% to 101% on the interpolation functions, the condition number given on 
Figure 5(b) [ remains close to the one computed with the Frobenius norm. Also, an important 
remark is that with K = 8 the computational effort for computing and is 

negligible compared to the one for and 



















K 

8 

16 

32 

64 

128 

256 

512 

Rescaled partial Hadamard 

0.4131 

0.3918 

0.3221 

0.1010 

0.0573 

0.0181 

0.0255 

Rescaled Rademacher (1) 

0.5518 

0.0973 

0.2031 

0.1046 

0.1224 

0.1111 

0.0596 

Rescaled Rademacher (2) 

1.0120 

0.6480 

0.1683 

0.1239 

0.0597 

0.0989 

0.0514 

Rescaled Rademacher (3) 

0.7193 

0.2014 

0.1241 

0.1051 

0.1235 

0.1369 

0.0519 

P-SRHT (1) 

0.4343 

0.2081 

0.2297 

0.0741 

0.0723 

0.0669 

0.0114 

P-SRHT (2) 

0.3624 

0.2753 

0.0931 

0.1285 

0.0622 

0.0619 

0.0249 

P-SRHT (3) 

0.8133 

0.4227 

0.1138 

0.0741 

0.0824 

0.0469 

0.0197 


Table 2: Relative error sup^ i|A(0 ~ II-^(0I|]R3- -^^(0 (resp. A(^)) are the in¬ 

terpolation functions associated to the Frobenius semi-norm projection (resp. the Frobenius 
norm projection) on Ym, with V either the rescaled partial Hadamard matrix, the random 
rescaled Rademacher matrix or the P-SRHT matrix (3 different samples for random matri¬ 
ces). 




(a) Interpolation functions. 


(b) Condition number of P 3 (^)T(^). 


Figure 5: Comparison between the Frobenius norm projection on Y 3 (black lines) and 
the Frobenius semi-norm projection on Kj, using for V either a sample of the rescaled 
Rademacher matrix (blue lines), the rescaled partial Hadamard matrix (red lines) or a sample 
of the P-SRHT matrix (green lines) with K = 8. 


5.1.3 Greedy selection of the interpolation points 

We now consider the greedy selection of the interpolation points presented in Section We 
start with an initial point = 0 and the next points are dehned by (38), where matrix V 
is a realization of the P-SRHT matrix with K = 128 columns. Pm{0 is the projection on 
Ym using the Frobenius semi-norm dehned by (|^. The hrst 3 steps of the algorithm are 
illustrated on Figure]^ We observe that at each iteration, the new interpolation point ^m+i 
is close to the point where the condition number of Pm(0^(0 i® maximal. Table presents 
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the maximal value over ^ G S of the residual, and of the condition number of Pm(0^(0- 
Both quantities are rapidly decreasing with m. This shows that this algorithm, initially 
designed to minimize ||(/ — Pm{0^i.0)^\\F, seems to be also efficient for the construction 
of preconditioners, in the sense that the condition number decreases rapidly. 



Figure 6: Greedy selection of the interpolation points: the hrst row is the residual ||(/ — 
(fhe blue points correspond to the maximum of the residual) with V a 
realization of the P-SRHT matrix with K = 128 columns, and the second row is the condition 
number of Pm(0^(0- 


iteration m 

0 

1 

2 

5 

10 

20 

30 


10001 

6501 

3037 

165,7 

51,6 

16,7 

7,3 

sup^||(/-P„(0A(0)F||f 

- 

300 

265 

80,5 

35,4 

10,0 

7,6 


Table 3: Convergence of the greedy algorithm: supremum over G S of the condition number 
(hrst row) and of the Frobenius semi-norm residual (second row). 


5.2 Multi-parameter-dependent equation 

We introduce a benchmark proposed within the OPUS project (see http://www.opus-project.fr). 
Two electronic components Vtjc (see Figure submitted to a cooling air how in the do¬ 
main VLait are hxed on a printed circuit board VtpcB- The temperature held dehned over 
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f2 = VLjc U fipsc U VLait C satisfies the advection-diffusion equation: 

-V■(/^(0Vu) + D(0^^■Vu = /. (44) 

The diffusion coefficient is equal to upcB on flpcB, f^Air on flAir and kjc on fljc- 
The right hand side / is equal to Q = 10® on Qjc and 0 elsewhere. The boundary conditions 
are n = 0 on T^, 62 ■ Vu = 0 on T^ (ei,e 2 are the canonical vectors of M^), and M|r, = M|r,. 
(periodic boundary condition). At the interface Tc = dfljc H dflpcB there is a thermal 
contact conductance, meaning that the temperature held u admits a jump over Tc which 
satishes 


Kic{ei ■ — K,pcB{ei ■ Vuinpcs) ~ ~ ^^pcb) 

The advection held v is given by v{x, y) = e 2 g{x), where g{x) = 0 if a; < epcB + e/c and 

2 a: — (e + ejc + “^gpcb) '' 


g{x) = 


2(e- e/c) 


ejc 


otherwise. We have 4 parameters: the width e := .^1 of the domain ^2Air, the thermal 
conductance parameter r := ^ 2 , the dihusion coefficient kjc '■= ^3 of the components and 
the amplitude of the advection held D := ,^ 4 . Since the domain = f2(e) depends on the 
parameter ,^1 G [emm,emax], we introduce a geometric transformation (x,r/) = 0 ^j(a:o,|/o) 
that maps a reference domain flo = ^(emax) to f 2 (^i): 


06 (^0,2/0) = 


with Co = epcB + ^ic- This method is described in [SB]: since the geometric transforma¬ 
tion 0^^ satishes the so-called Affine Geometry Precondition, the operator of equation (44) 
formulated on the reference domain admits an affine representation. 

For the spatial discretization we use a hnite element approximation with n = 2.8 ■ 10'^ 
degrees of freedom (piecewise linear approximation). We rely on a Galerkin method with 
SUPG stabilization (see i)- S is a set of 10 ^ independent samples drawn according the 
loguniform probability laws of the parameters given on Figure 


(\ 

Xo 

if Xo < eo 


t 

eo + (^0 — eo)-^^— 

^ ^ ^ ^ emax-eic 

otherwise. 


V 

Vo 


/ 


5.2.1 Preconditioner for the projection on a given reduced space 

We consider here a POD basis Xr of dimension r = 50 computed with 100 snapshots of 
the solution (a basis of is obtained by the hrst 50 dominant singular vectors of a matrix 
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Geometry (m) 

II 

~ logi/(2.5 X 10“y5 X 10-q 

epCB 

2 X 10“^ 

eic 

2 X 10“^ 

hic 

2 X 10“^ 

hi 

3 X 10“^ 

^2 

2 X 10“^ 


4 X 10“^ 

Thermal conductance (T^m '‘K 

r = i2 

~ logU{l X IQ-yi X 10^) 

Diffusion coefficients [Wm 

KPCB 

2 X 10“-^ 

KAir 

3 X 10“^ 

II 

0 

~ \ogU{2 X 10“51.5 X lOh 

Advection field {Wm ‘‘K 

0=^4 

~ logU{5 X lO*'*,! X 10“^) 



Figure 7: Geometry and parameters of the benchmark OPUS. 


of 100 random snapshots of m( 0)- Then we compute the Petrov-Galerkin projection as 

The efficiency of the preconditioner can be measured with the 
(.^): the associated quasi-optimality constant (1 — should be as 


presented in Section |3.1 
quantity h 


close to one as possible (see equation (33)). We introduce the quantile Qp of probability p 
associated to the quasi-optimality constant (1 — defined as the smallest value 

gp > 1 satisfying 

P({^ e S : (1 - < qp}) > P, 


where P(74) = ^A/^'E for A<zE. Table shows the evolution of the quantile with respect 
to the number of interpolation points for the preconditioner. Here the goal is to compare 
the different strategies for the selection of the interpolation points: 


(a) the greedy selection (42) based on the quantity hr,m(0) 


(b) the greedy selection (38) based on the Frobenius semi-norm residual, with V a P-SRHT 
matrix with K = 256 columns, and 


(c) a random Latin Hypercube sample (LHS). 

The projection on (or Y^) is then defined with the Frobenius semi-norm using for V a 
P-SRHT matrix with K = 330 columns. 

When considering a small number of interpolation points m < 3, the projection on Y^ 
provides lower quantiles for the quasi-optimality constant compared to the projection on Wi- 
The positivity constraint is useful for small m. But for high values of m (see m = 15) the 
positivity constraint is no longer necessary and the projection on Ym provides lower quantiles. 

Goncerning the choice of the interpolation points, the strategy (a) shows the faster decay 
of the quantiles Qp, especially for p = 50%. The strategy (b) shows also good results, but the 
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quantile Qp for p = 100% are still high compared to (a). These results show the benefits of 
the greedy selection based on the quasi-optimality constant. Finally the strategy (c) shows 
bad results (high values of the quantiles), especially for small m. 







Projection 

on 







Greedy selection based on 


(c) Latin Hypercube 



(a) %.m(G 

(b) Frob. residual 


sampling 



50% 

90% 

100% 

50% 

90% 

100% 

50% 

90% 

100% 

m 

= 0 

21.3 

64.1 

94.1 

21.3 

64.1 

94.1 

21.3 

64.1 

94.1 

m 

= 1 

18.3 

74.1 

286.7 

10.2 

36.1 

161.6 

18.3 

104.1 

231.8 

m 

= 2 

11.9 

22.6 

42.1 

9.8 

53.3 

374.0 

11.5 

113.0 

533.9 

m 

= 3 

11.1 

49.2 

200.4 

7.8 

31.2 

60.2 

18.3 

138.7 

738.5 

m 

= 5 

5.2 

10.8 

18.4 

6.8 

18.6 

24.5 

8.7 

121.1 

651.4 

m 

= 10 

3.1 

9.0 

13.2 

5.3 

22.3 

62.1 

4.0 

21.6 

345.7 

m 

= 15 

2.2 

6.3 

10.4 

3.5 

6.5 

11.5 

2.7 

7.8 

48.6 







Projection 

onH+ 







Greedy selection based on 


(c) Latin Hypercube 



(a) 5r,m{.0 

(b) Frob. residual 


sampling 



50% 

90% 

100% 

50% 

90% 

100% 

50% 

90% 

100% 

m 

= 0 

21.3 

64.1 

94.1 

21.3 

64.1 

94.1 

21.3 

64.1 

94.1 

m 

= 1 

18.3 

74.1 

286.7 

10.2 

36.1 

161.6 

18.3 

104.1 

231.8 

m 

= 2 

11.9 

22.6 

42.1 

8.9 

35.5 

78.6 

10.4 

41.5 

112.5 

m 

= 3 

9.7 

24.4 

48.0 

7.9 

27.7 

57.9 

12.1 

48.8 

114.1 

m 

= 5 

6.4 

15.0 

25.5 

6.9 

26.8 

65.1 

5.7 

11.6 

17.5 

m 

= 10 

4.6 

9.5 

16.8 

7.3 

18.9 

38.0 

4.3 

10.0 

18.5 

m 

= 15 

4.3 

7.1 

11.2 

6.4 

10.1 

18.0 

4.2 

9.0 

19.3 


Table 4: Quantiles Qp of the quasi-optimality constant associated to the Petrov-Galerkin 
projection on the POD subspace Xr for p = 50%, 90% and 100%. The row m = 0 corresponds 
to Pb(0 = is the standard Galerkin projection. 


5.2.2 Preconditioner for Reduced Basis method 

We now consider the preconditioned Reduced Basis method for the construction of the ap¬ 
proximation space Xr, as presented in Section [3^ Figuresand [T^ show the convergence of 
the error with respect to the rank r of Ur{0 different constructions of the preconditioner 
PmiO- Two measures of the error are given: sup^g= ||'*^(0 ~ w(0llx/||w(0llx, and the quan¬ 
tile of probability 0.97 for i|M(0 ~'Wr(0llx:/||w(0llx- The curve “Ideal greedy” corresponds 
to the algorithm defined by ( p4| ) which provides a reference for the ideally conditioned algo¬ 
rithm, i.e. with KmiO — 1- Figure shows the corresponding hrst interpolation points for 
the solution. 
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The greedy selection of the interpolation points based on (38) (see Figure]^ allows to 
almost recover the convergence curve of the ideal greedy algorithm when using the projection 
on Ym with m = 15. For the strategy of re-using the operators factorizations, the approxi¬ 
mation is rather bad for r = m < 10 meaning that the space Fj. (or ) is not really adapted 
for the construction of a good preconditioner over the whole parametric domain. However, 
for higher values of r, the preconditioner is getting better and better. For r > 20, we al¬ 
most reach the convergence of the ideal greedy algorithm. We conclude that this strategy of 
re-using the operator factorization, which has a computational cost comparable to the stan¬ 
dard non preconditioned Reduced Basis greedy algorithm, allows obtaining asymptotically 
the performance of the ideal greedy algorithm. Note that the positivity constraint yields a 
better preconditioner for small values of r but is no longer necessary for large r. 





(a) (b) (c) (d) (e) (f) 



(a) r = 1 

(b) r = 2 

(c) r = 3 

(d) r = 4 

(e) r = 5 

(f) r = 6 

e 

7.1 • 10-3 

6.1 • 10-3 

5.0-10-3 

5.0-10-2 

5.1 - 10-3 

4.4 - 10-2 

r 

9.9- IQi 

2.3- 10° 

4.1 • 10-3 

2.8- 10° 

8.6 - 10-3 

4.8 - 103 

KIC 

1.1 • 10^ 

2.1 • 10-3 

3.2 • 103 

6.0- 10° 

1.1 -102 

3.2 - 10-3 

D 

1.7- 10-3 

7.4 • 10-4 

6.2 • 10-4 

9.9-10-3 

8.6 -10-3 

6.0 - 10-4 


Figure 8: First six interpolation points of the ideal reduced basis method and corresponding 
reduced basis functions. 

Let us hnally consider the effectivity index 

h.(0 = iin(0(^(eK(e) -&(0)iu/ii«(o - w(e)iix, 

which evaluates the quality of the preconditioned residual norm for error estimation. We 
introduce the conhdence interval Ir{p) dehned as the smallest interval which satishes 

P({^ e S : rir{i) e Ir{,p)}) >P- 

On Figure pT] we see that the conhdence intervals are shrinking around 1 when r increases, 
meaning that the preconditioned residual norm becomes a better and better error estimator 
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(a) Projection on (b) Projection on 




(c) Projection on Ym 


(d) Projection on Y^ 


Figure 9: Convergence of the preconditioned reduced basis method using the greedy selection 
of interpolation points for the preconditioner. Supremum over S (top) and quantile of 
probability 97% (bottom) of the relative error ||m( 0 ~ '^r(Ollx/||t^('C)llx with respect to 
r. Comparison of preconditioned reduced basis algorithms with ideal and standard greedy 
algorithms. 


when r increases. Again, the positivity constraint is needed for small values of r, but 
we obtain a better error estimation without imposing this constraint for r > 20. On the 
contrary, the standard residual norm leads to effectivity indices that spread from 10“^ to 
10 ^ with no improvement as r increases, meaning that we can have a factor 10^ between the 
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Figure 10: Preconditioned Reduced basis methods with re-use of operators. Supremum over 
S (left) and quantile of probability 97% (right) of the relative error \\u{^) —Mr(0llv/||w(0llv 
with respect to r. Comparison of preconditioned reduced basis algorithms with ideal and 
standard greedy algorithms. 


error estimator ||%(^)Mr(0 ~ ^(OiU' and the true error ||Mr(0 ~ '^(Oliv¬ 



ia) Re-use, proj. on Ym- (b) Re-use, proj. on Y^. (c) Dual residual norm. 


Figure 11: Conhdence intervals of the effectivity index during the iterations of the Reduced 
Basis greedy construction. Comparison between preconditioned algorithms with re-use of 
operators factorizations (a,b) and the non preconditioned greedy algorithm (c). 
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6 Conclusion 


We have proposed a method for the interpolation of the inverse of a parameter-dependent 
matrix. The interpolation is dehned by the projection of the identity in the sense of the 
Frobenius norm. Approximations of the Frobenius norm have been introduced to make com¬ 
putationally feasible the projection in the case of large matrices. Then, we have proposed 
preconditioned versions of projection-based model reduction methods. The preconditioner 
can be used to dehne Petrov-Galerkin projections on a given approximation space with bet¬ 
ter quasi-optimality constants by introducing a parameter-dependent test space depending 
on the preconditioner. Also, the preconditioner can be used to improve residual-based er¬ 
ror estimates that are used for assessing the quality of a given approximation, which is 
required in any adaptive approximation strategy. Different strategies have been proposed 
for the selection of interpolation points depending on the objective: (i) the construction of 
an optimal approximation of the inverse operator for preconditioning iterative solvers or for 
improving error estimators based on preconditioned residuals, (ii) the improvement of the 
quality of Petrov-Galerkin projections of the solution of a parameter-dependent equation 
on a given reduced approximation space, or (iii) the re-use of operators factorizations when 
solving a parameter-dependent equation with a sample-based approach. The performance 
of the obtained parameter-dependent preconditioners has been illustrated in the context of 
projection-based model reduction techniques such as the Proper Orthogonal Decomposition 
and the Reduced Basis method. 

The proposed preconditioner has been used to dehne Petrov-Galerkin projections with 
better stability constants. For the solution of PDFs, the Petrov-Galerkin projection has 
been dehned at the discrete (algebraic) level for obtaining a better approximation (in a 
reduced space) of the hnite element Galerkin approximation of the PDF. Therefore, for 
convection-dominated problems, the proposed approach does not avoid using stabilized hnite 
element formulations. Similar observations can be found in [3l]. However, a Petrov-Galerkin 
method could be dehned at the continuous level with a preconditioner being the interpolation 
of inverse partial diherential operators. In this continuous framework, the preconditioner 
would improve the stability constant for the hnite element Galerkin projection and may 
avoid the use of stabilized hnite element formulations. Such Petrov-Galerkin methods have 
been proposed in [niDS] for convection-dominated problems (as an alternative to standard 
stabilization methods), which can be interpreted as an implicit preconditioning method 
dehned at the continuous level. 

In the present paper, the parameter-dependent preconditioner is obtained by a projection 
onto the space generated by snapshots of the inverse operator. When the storage of many 
inverse operators (even as implicit matrices) is not feasible, a parameter-dependent precon- 
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ditioner could be obtained by a projection into the linear span of preconditioners, such as 
incomplete factorizations, sparse approximate inverses, H-matrices or other preconditoners 
that are readily available for a considered application. Also, we have restricted the presenta¬ 
tion to the case of real matrices but the methodology can be naturally extended to the case 
of complex matrices. 
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